/*
 * KineticsSolver.cpp
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#include "KineticSolverEuler.h"

// Eigen includes
#include <Eigen/LU>

// GeoReact includes
#include "Utils/Algorithms.h"
#include "Utils/Assert.h"

KineticSolverEuler::KineticSolverEuler()
{}

KineticSolverEuler::KineticSolverEuler(const Multiphase& multiphase, 
                                       const ReactionSystem& eReactions, 
                                       const ReactionSystem& kReactions,
                                       const vector<string>& jSpecies, 
                                       const vector<string>& eSpecies, 
                                       const vector<string>& kSpecies)
{
	Initialiaze(multiphase, eReactions, kReactions, jSpecies, eSpecies, kSpecies);
}

KineticSolverEuler::~KineticSolverEuler()
{}

void KineticSolverEuler::Initialiaze(const Multiphase& multiphase, 
                                     const ReactionSystem& eReactions, 
                                     const ReactionSystem& kReactions,
                                     const vector<string>& jSpecies, 
                                     const vector<string>& eSpecies, 
                                     const vector<string>& kSpecies)
{
	// Set data members
	this->multiphase = multiphase;
	this->eReactions = eReactions;
	this->kReactions = kReactions;
	
	// Set the number of primary (Nj), equilibrium (Ne) and kinetic (Nk) species
	Nj = jSpecies.size();
	Ne = eSpecies.size();
	Nk = kSpecies.size();
	
	// Set the indexes of the primary (jIndexes), equilibrium (eIndexes) and kinetic (kIndexes) species
	jIndexes = multiphase[jSpecies];
	eIndexes = multiphase[eSpecies];
	kIndexes = multiphase[kSpecies];
	
	// Assemble the stoichiometric matrices (ukj) and (ukk)
	ukj = kReactions.AssembleStoichiometricMatrix(jSpecies);
	ukk = kReactions.AssembleStoichiometricMatrix(kSpecies);
	
	// Assemble the stoichiometric matrices (uej) and (uee)
	MatrixXd uej = eReactions.AssembleStoichiometricMatrix(jSpecies);
	MatrixXd uee = eReactions.AssembleStoichiometricMatrix(eSpecies);
	
	// Compute (uee_inv), which is the inverse of the stoichiometric matrix (uee)
	FullPivLU<MatrixXd> lu(uee);
	
	MatrixXd uee_inv = lu.inverse();
	
	// Compute the canonical stoichiometric matrix of the equilibrium-controlled reactions
	vej = - uee_inv * uej;
	
	// Initialise the auxiliar Eigen vectors with zero entries
	nj = VectorXd::Zero(Nj);
	ne = VectorXd::Zero(Ne);
	nk = VectorXd::Zero(Nk);
	uj = VectorXd::Zero(Nj); 
	
	// Initialiaze the equilibrium solver
	eSolver.Initialiaze(multiphase, eReactions, jSpecies, eSpecies);
}

bool KineticSolverEuler::Iterate(double& t, double tEnd, double& dt, double T, double P, VectorXd& n)
{
	// Transfer the nj, ne, and nk part of the vector (n) to the vectors (nj), (ne), and (nk) respectively
	for(unsigned j = 0; j < Nj; ++j) nj[j] = n[jIndexes[j]];
	for(unsigned e = 0; e < Ne; ++e) ne[e] = n[eIndexes[e]];
	for(unsigned k = 0; k < Nk; ++k) nk[k] = n[kIndexes[k]];
	
	// Calculate the total composition vector (uj), where uj = nj + [vej]T * ne
	uj = nj + vej.transpose() * ne;
	
	// Calculate the activities of the species, now that the composition vector (n) has been determined
	const VectorXd a = multiphase.Activities(T, P, n);
	
	// Calculate the mass rates of the kinetic-controlled reactions with temperature (T), pressure (P), composition vector (n), and activity vector (a)
	const VectorXd rk = kReactions.KineticRates(T, P, n, a);
	
	// Update (uj) and (nk) through the ordinary differential equations
	uj += dt * ukj.transpose() * rk;
	nk += dt * ukk.transpose() * rk;
	
	// Update the nk part of the vector (n)
	for(unsigned k = 0; k < Nk; ++k) n[kIndexes[k]] = nk[k];
	
	// Solve the equilibrium problem with total composition of the primary species (uj)
	auto diagnostic = eSolver.Solve(uj, T, P, n);
	
	Assert(diagnostic.converged, "Method did not converged.");
	
	t += dt;
	
	return diagnostic.converged;
}

